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Preface 


The  study  described  herein  was  conducted  during  1999  for  the  U.S.  Army 
Engineer  District,  Jacksonville,  and  the  St.  Johns  River  Water  Management 
District,  Palatka,  FL,  by  personnel  of  the  Coastal  and  Hydraulics  Laboratory 
(CHL)  of  the  U.S.  Army  Engineer  Research  and  Development  Center  (ERDC), 
Vicksburg,  MS.  General  supervision  was  provided  by  Dr.  James  R.  Houston, 
former  Director  of  CHL;  Mr.  Thomas  W.  Richardson,  Acting  Director,  CHL;  and 
Dr.  Robert  T.  McAdory,  Chief  of  Tidal  Hydraulics  Branch  (HE-T),  CHL. 

Messrs.  Jose  A.  Sanchez,  of  CHL,  and  Robert  A.  Evans,  Jr.  of  the  Jacksonville 
District  conducted  the  study  and  prepared  this  report. 

At  the  time  of  publication  of  this  report,  Dr.  James  R.  Houston  was  Director  of 
ERDC  and  COL  Robin  R.  Cababa,  EN,  was  Commander. 

This  report  should  be  cited  as  follows: 

Sanchez,  Jose  A.,  and  Evans,  Robert  A.  Jr.  (2000).  “The  Lower  St.  Johns 
River  and  Estuary  Numerical  Model  Study,”  Technical  Report  ERDC/CHL 
TR-00-10,  U.S.  Army  Engineer  Research  and  Development  Center,  Vicksburg, 
MS. 


The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or 
promotional  purposes.  Citation  of  trade  names  does  not  constitute  an  official 
endorsement  or  approval  of  the  use  of  such  commercial  products. 


1  Introduction 


The  U.S.  Army  Engineer  Research  and  Development  Center  (ERDC)  was 
asked  to  develop  a  multidimensional  hydrodynamic  and  salinity  model  for  the 
Lower  St.  Johns  River  and  Estuary  (LSJRE),  Florida.  The  developed  model  was 
intended  as  a  general  planning  tool  to  gain  insight  into  the  hydrodynamics  and 
transport  mechanisms  of  the  system. 

The  LSJRE  has  some  salinity  stratification  in  some  of  its  lower  reaches  and  is 
well  mixed  in  the  upper  reaches.  To  adequately  duplicate  this  stratification  and  the 
associated  hydrodynamics,  a  multidimensional  model  is  needed.  In  a  previous 
effort,'  the  modeling  approach  was  recommended  based  on  expected  model  needs. 
This  report  describes  the  implementation  of  those  recommendations. 


Project  Objectives 


The  objective  of  this  study  was  to  develop  a  verified  multidimensional  model 
of  hydrodynamics  and  salinity  transport  in  the  Lower  St.  Johns  River  and 
Estuary.  The  work  was  accomplished  in  the  following  steps: 

a.  Develop  a  working,  yet  unverified,  multidimensional  model  of 
hydrodynamics  and  salinity  transport  in  the  LSJRE. 

b.  Verify  the  multidimensional  model  of  hydrodynamics  and  salinity 
transport  against  field  measurements  of  tidal  elevation,  velocity,  and 
salinity  concentrations. 

Task  1  was  the  model  development  phase  of  the  project.  Model  development 
encompasses  defining  questions  to  be  posed,  creating  the  numerical  mesh,  testing 
the  mesh  using  synthetic  boundary  conditions,  refining  and  reworking  of  the  mesh 
to  improve  model  accuracy  and  stability,  and  assigning  initial  model  parameters. 
To  improve  performance,  modifications  of  the  numerical  code  were  included. 

Task  2  was  the  verification  phase  of  the  project.  Verification  includes: 


1  Roig,  Lisa  C.  (1994).  “Review  and  evaluation  of  hydrodynamic  modeling  for  the  Lower  St.  Johns  River 
Estuary,”  Technical  Report  HL-94-15,  U.S.  Army  Engineer  Waterways  Experiment  Station,  Vicksburg,  MS. 
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a.  the  development  of  synoptic,  comprehensive  field  (prototype)  data 
sets  for  boundary  conditions  and  model  adjustment, 

b.  model  spin-up  to  achieve  realistic  initial  conditions, 

c.  performance  of  model  runs, 

d.  comparison  of  simulation  results  against  prototype  data, 

e.  analysis  of  simulation  results, 

f.  adjustments  to  improve  model  performance, 

g.  comparison  of  the  adjusted  model  against  prototype  data, 

h.  repeat  steps  e  through  g  until  satisfactory  results  are  obtained. 

Lower  St.  Johns  River  and  Estuary  Description 


The  Lower  St.  Johns  River  and  Estuary  extends  from  the  mouth 
approximately  161  km  (100  miles),  to  just  north  of  Welaka,  FL.  The  mean  range 
of  the  tide  is  1.67  m  (5.49  ft)  at  the  mouth,  and  0.1 1  m  (0.35  ft)  at  Welaka. 

Between  these  two  points,  the  mean  tide  range  varies  from  0.22  m  (0.71  ft)  at 
Julington  Creek  to  0.33  m  (1.09  ft)  at  Palatka.  These  differences  are  due  mainly 
to  considerable  variations  of  the  channel  width. 

The  salinity  distribution  varies  depending  on  the  distance  from  the  ocean.  At 
the  ocean,  the  salinity  is  between  30  -  35  ppt;  at  Buffalo  Bluff,  the  salinity  is  less 
than  0.5  ppt.  The  vertical  salinity  distribution  varies  between  a  well  mixed 
structure  and  a  slightly  stratified  structure.  Based  on  the  data  collected  June  1995 
-  February  1996,  the  stretch  of  the  river  above  the  Shands  Bridge  is  well  mixed; 
below  the  Buckman  Bridge  there  are  periods  when  the  salinity  is  stratified.  In 
between  these  locations,  it  is  assumed  that  the  salinity  varies  between  stratified 
and  well  mixed. 

The  1994  bathymetric  survey  covered  only  the  main  stem  of  the  river. 
Bathymetric  data  for  the  tributaries  were  taken  from  the  National  Oceanic  and 
Atmospheric  Administration  (NOAA)  charts,  which  is  referenced  to  mean  lower 
low  water  (mllw).  The  nodes  of  the  initial  finite  element  mesh  for  the  tributaries 
were  digitized  from  these  maps.  The  Surface  Water  Modeling  System  (SMS) 
Graphical  Users  Interface  (GUI)  was  then  used  to  unite  the  tributary  meshes  with 
the  main  stem  mesh.  Discrepancies  among  the  datums  of  the  bathymetric  data 
from  the  various  sources  were  carefully  corrected.  The  bathymetry  of  the  model 
was  referenced  to  the  National  Geodetic  Vertical  Datum  (NGVD).  The  difference 
between  mllw  and  NGVD  was  obtained  from  published  data-  (it  must  be  noted 
that  the  published  data  were  preliminary). 

The  character  of  the  St.  Johns  River  itself  influenced  many  of  the  mesh 
decisions.  For  example,  a  large  expanse  of  tidal  marshes  exists  to  the  north  of  the 
river  near  its  mouth.  These  tidal  marshes  are  known  to  store  and  release  tidal 
waters  entering  and  exiting  the  marshes  through  several  inlets.  The  volume  of  this 
tidal  exchange  is  unknown  and  difficult  to  measure. 


2  National  Oceanic  and  Atmospheric  Administration  (1993).  “Final  Report,  St.  Johns  River,  Florida,  Tide  control 
and  tidal  characteristics''. 
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In  a  two-dimensional,  channel  deepening  study  of  the  Lower  St.  Johns  River 
conducted  by  ERDC  for  the  U.  S.  Army  Engineer  District,  Jacksonville  during 
1994  and  1995,  it  was  noted  that  the  simulation  model  had  difficulty  matching 
velocity  measurements  in  the  vicinity  of  the  entrance  to  the  Intracoastal  Waterway. 
It  was  thought  that  the  major  factors  affecting  these  currents  were  the  tidal 
volume  stored  in  the  adjacent  marshes  and  the  modulation  of  the  tidal  signal 
resulting  from  this  storage.  A  numerical  experiment  was  carried  out  during  this 
earlier  study  bearing  out  these  assumptions. 

For  this  modeling  effort,  the  magnitude  of  the  water  discharge  interchange 
between  the  main  stem  of  the  river  and  the  tidal  marshes  through  two  tributaries, 
Browns  Creek  and  Clapboard  Creek,  was  estimated  using  discharge  data  collected 
during  1994.  The  same  was  done  for  the  northern  and  southern  sections  of  the 
Intracoastal  waterway. 

A  second  characteristic  of  the  St.  Johns  River  influencing  the  level  of  mesh 
resolution  is  the  number  of  tributaries  and  inlets.  For  reasons  of  model  accuracy, 
it  is  necessary  to  increase  the  resolution  around  bends,  comers,  peninsulas, 
islands,  channel  expansions  and  channel  contractions  to  maintain  smooth  mesh 
transitions.  Geometric  features  that  have  a  significant  influence  on  the  overall 
circulation  of  the  estuary  should  be  adequately  resolved  to  represent  the  flow 
patterns  caused  by  the  feature.  If  the  feature  does  not  have  a  significant  influence 
on  the  overall  circulation,  then  it  is  better  to  ignore  the  feature  and  preserve  the 
smoothness  of  the  mesh.  Irregular  regions  of  the  mesh  can  cause  spurious 
oscillations  of  the  numerical  model  solution  and  degrade  solution  accuracy.  The 
sinuosity  of  the  main  stem,  the  numerous  branch  points  and  inlets,  and  the 
numerous  large  tributaries  contribute  to  the  complex  geometry  of  this  system. 

Thus,  the  level  of  resolution  necessarily  increases  with  the  number  of  features 
included.  Unfortunately,  as  the  level  of  resolution  is  increased,  the  computer  time 
required  for  processing  increases. 

It  was  determined  that  the  mesh  must  extend  into  the  main  tributaries  for 
several  reasons.  First,  inflow  boundary  conditions  must  be  specified  at  each  of  the 
major  tributaries  to  properly  account  for  the  mixing  of  fresh  and  brackish  waters 
having  more  saline  water  masses.  Boundary  locations  were  moved  upstream  to 
points  along  the  tributaries  where  inflow  salinities  were  thought  to  be  known.  The 
St.  Johns  River  Water  Management  District  (SJRWMD)  specified  the  locations  of 
the  boundary  data  points.  Second,  the  circulation  patterns  of  the  main  stem  in  the 
vicinities  of  the  tributaries  is  generally  unknown.  The  solution  domain  must  be 
contiguous  through  the  zone  of  influence  of  the  main  stem  to  properly  calculate 
the  tidal  circulation. 
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2  Numerical  Hydrodynamic/Salinity 
Model 


RMA10-WES 


The  generalized  numerical  hydrodynamic  model  used  in  this  project  is  known 
as  RMA10-WES.  The  model  was  originally  developed  under  contract  to 
Waterways  Experiment  Station,  now  known  as  ERDC,  by  Dr.  Ian  King  of 
Resource  Management  Associates  of  Suisun  City,  CA.  The  version  known  as 
RMA10-WES  has  been  modified  significantly  by  ERDC  personnel  and  is 
described  in  Appendix  A.  A  summary  of  technical  specifications  of  the  model  is 
given  below: 

a.  Galerkin  finite  element  formulation, 

b.  Unstructured  mesh  composed  of  bricks,  tetrahedral  elements,  prism, 
and  pyramids  in  3-D,  triangles  and  quadrilaterals  in  2-D  (horizontal 
averaging  or  vertical  averaging),  and  line  elements  in  1-D; 

c.  Mixed  interpolation  scheme  for  pressures  and  velocities  (linear  and 
quadratic); 

d.  Z  -  based  vertical  coordinate  transform; 

e.  Governing  equations  of  the  nonlinear  Reynolds  form  of  the  Navier- 
Stokes  equations,  including  baroclinic  forcing,  wind  forcing,  tidal 
forcing,  Coriolis  effect,  and  bed  friction; 

f.  Three-dimensional,  two-dimensional  vertically  averaged,  two- 
dimensional  laterally  averaged,  and  one-dimensional  approximations 
within  the  same  numerical  mesh; 

g.  Time  integration  via  an  implicit  Crank-Nicholson  finite  difference 
operator; 

h.  Eddy  viscosity  formulation  for  horizontal  turbulence  closure; 

i.  Mellor-Yamada  Level  2  vertical  turbulence  closure  model; 

j.  Non-linear  acceleration  and  friction  terms  of  the  governing  equations; 

k.  Newton-Raphson  iteration  for  the  solution  of  the  non-linear  system  of 
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equations; 

l.  Frontal  solver  for  the  solution  of  the  matrix  equations;  and 

m.  Two  algorithms  available  for  wetting  and  drying  of  2-D  vertically 
averaged  elements. 


Numerical  Model  Mesh 


The  numerical  model  mesh  boundary  of  the  Lower  St.  Johns  River  is  shown 
in  Figures  la  to  le.  The  elements  of  the  mesh  were  divided  into  57  categories, 
called  material  types.  In  general,  these  material  types  were  defined  by  size, 
location,  and  average  depth.  Table  1  lists  the  material  types  and  the  associated 
hydraulic  and  transport  parameters.  Table  1  also  lists  the  number  of  layers 
specified  for  each  material  type  when  the  model  is  used  in  3-D  mode. 

The  ocean  area  is  classified  with  material  types  2,  15,  16,  and  17.  When  run 
in  3-D  mode,  all  of  these,  with  the  exception  of  Material  15,  are  three-dimensional 
elements.  Material  15  is  specified  as  two-dimensional,  because  it  is  on  the  ocean 
boundary.  Specifying  the  ocean  tidal  boundary  along  two-dimensional  elements 
aids  in  model  stability.  In  addition,  eddy  viscosity  and  Manning's  n  values  in  the 
ocean  are  artificially  high  to  control  instability  near  the  boundaries.  These  values 
are  numerically,  rather  than  physically  based,  but  do  not  compromise  results  in  the 
interior  part  of  the  model.  Additionally,  the  region  just  outside  of  the  jetties  has 
been  schematized  to  insure  a  smooth  boundary.  The  reason  for  not  including  the 
jetties  precisely  is  to  reduce  the  front  width  of  the  model,  and  thereby  the 
execution  time. 

The  numerical  mesh  is  defined  by  5696  1-D/2-D  elements  and  16457  l-D/2- 
D  nodes.  Extending  the  model  to  three-dimensions  increases  the  total  number  to 
15,780  elements  and  37,187  nodes.  The  three-dimensional  region  of  the  model 
only  includes  the  main  stem  up  to  Shands  Bridge  and  part  of  the  ocean.  In  areas 
where  the  stratification  is  significant,  the  three-dimensional  area  is  defined  by  two, 
three,  or  four  layers;  in  areas  where  the  stratification  is  less  significant,  there  is 
only  one  layer;  and,  in  areas  where  there  is  no  stratification,  the  elements  are  two- 
dimensional.  For  finite-element,  quadratic  elements  of  the  type  used  in  RMA10- 
WES,  the  number  of  vertical  computation  nodes  is  greater  than  the  number  of 
layers.  The  number  of  vertical  computation  nodes  in  an  area  is  found  using  the 
following  formula: 


Ny  =  2  NL  +  1 


where 

Ny  =  number  of  vertical  computation  nodes, 

Nl  =  number  of  layers  specified  in  the  area. 

Some  of  the  minor  tributaries  and  marsh  areas  are  represented  by  one- 
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dimensional  elements.  This  approach  minimizes  the  computer  simulation  time 
requirements. 


Boundary  Conditions 


Figure  la  shows  the  location  of  nine  inflow  data  stations  used  as  model 
boundary  conditions.  The  SJRWMD  provided  average  daily  discharge  values  for 
Big  Davis  Creek,  Black  Creek,  Buffalo  Bluff,  Cedar  River,  Deep  Creek,  Dunns 
Creek,  Ortega  River,  Rice  Creek,  and  Six  Mile  Creek.  The  data  covered  the 
period  1  April  1996  to  31  January  1997.  Recording  of  discharge  data  at  Cedar 
River  stopped  on  30  September  1996.  Discharge  records  of  Deep  Creek,  Dunns 
Creek,  Ortega  River,  and  Six  Mile  Creek  showed  several  gaps  in  the  data  starting 
during  the  first  days  of  October  1996.  Gaps  were  also  found  in  the  salinity  and 
water  surface  elevation  data  provided.  Based  on  the  system  data  on  hand  and  in 
consultation  with  the  SJRWMD,  120  days  were  selected,  starting  on  29  May  1996 
as  the  verification  period  to  avoid  major  gaps  in  the  field  data  driving  the  model. 

The  cumulative  total  stream  inflow  to  the  model  is  shown  in  Figure  2.  The  period 
represented  begins  on  29  April  1996  and  ends  on  26  September  1996.  Thirty 
days  were  added  to  the  beginning  of  the  verification  period  as  spin-up  time  to 
eliminate  any  residual  effects  of  the  initial  conditions  assigned  during  the  start  of 
the  simulation.  The  actual  flow  rates  are  shown  in  Figures  3-11.  The  primary 
inflow  sources  are  Dunns  Creek  and  Buffalo  Bluff.  The  discharge  at  Buffalo 
Bluff  contributes  more  inflow  than  all  other  streams  combined. 

Still,  flow  data  for  Cedar  River  and  Six  Mile  Creek  were  incomplete.  No 
flow  measurements  were  collected  for  Cedar  River  between  15  June  and  9  July. 

Forty-four  days  worth  of  data,  starting  on  29  May,  were  missing  from  the  Six  Mile 
Creek  data.  Gaps  were  filled  with  flow  values  derived  from  the  rest  of  the 
tributaries.  Daily  factors  were  calculated  for  each  tributary  based  on  the  formula: 

daily  factor  =  (daily  average  flow)  /  (average  flow  for  the 

verification  period,  excluding  days  when  data  was 
missing  at  Cedar  River  or  Six  Mile  Creek). 

Average  daily  factors  were  calculated  using  all  tributaries  and  were 
then  applied  to  estimate  daily  flows  at  Cedar  River  and  Six  Mile  Creek. 

The  tidal  head  (water  level)  boundary  is  shown  in  Figure  la.  The  tidal 
elevations  are  applied  along  the  northern  and  southern  edges  of  the  ocean 
boundary  of  the  model,  corresponding  to  the  area  in  front  of  Femandina  Beach 
and  Jacksonville  Beach  respectively.  The  elevations  specified  are  equal  and  in 
phase.  The  elevations  along  the  boundary  are  representative  of  the  tide  at  U.S. 

Coast  Guard  Base  at  Mayport.  Real  tide  data  measured  every  hour  at  Mayport 
was  provided  by  the  SJRWMD  and  was  applied  using  half  hour  time  steps  at  the 
model  boundary,  after  increasing  the  tide  range  by  1 8  percent.  This  amplification 
of  the  tide  range  was  necessary  to  account  for  the  energy  loss  between  the 
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boundary  and  Mayport.  A  wave  travel  time  of  1 5  minutes  between  the  tide 
boundary  and  Mayport  was  calculated.  The  travel  time  was  subtracted  from  the 
model  time  after  extracting  results  for  the  sake  of  compatibility  with  the  prototype. 
The  salinity  at  the  ocean  was  set  at  a  constant  value  of  35  ppt. 

Originally,  predicted  tide  signals  for  the  tide  boundary  conditions  at 
Jacksonville  Beach  and  Femandina  Beach  were  to  be  used.  This  approach  would 
have  yielded  variation  between  the  two  tidal  boundaries  to  more  accurately  drive 
the  currents  in  the  ocean  area  of  the  mesh.  The  difference  between  predicted  and 
measured  tide  at  Mayport  was  calculated  to  account  for  atmospheric  events  not 
included  in  the  predicted  tide  signal.  This  difference  was  then  applied  to  the 
predicted  tide  signals  used  as  boundary  condition  at  Jacksonville  Beach  and 
Femandina  Beach.  A  long-term  simulation  (five  months)  showed  that  the  modeled 
tide  signal  diverged  from  the  prototype.  This  fact  repeated  itself  at  all  tide  stations 
used  for  verification.  The  problem  was  caused  by  a  defect  in  the  predicted  tide 
data.  After  plotting  the  daily  average  tide  during  the  five  month  period,  a 
downward  curve  was  observed.  As  a  consequence,  this  approach  was  discarded. 

The  ocean  region  of  the  mesh  is  the  area  of  the  model  outside  the  jetties.  It 
extends  approximately  18  km  (1 1  miles)  east  beyond  the  end  of  the  jetties  and 
approximately  14.5  km  (9  miles)  north  and  south  of  the  jetties.  These  dimensions 
were  required  to  minimize  occurrences  of  fresher  water  from  the  river  reaching 
the  ocean  boundary,  which  could  have  caused  numerical  instabilities. 

In  addition  to  stream  flow  and  tidal  boundaries,  wind,  rainfall,  and 
evaporation  were  considered.  Historically,  the  period  of  May- August  has  more 
rainfall  than  other  seasons.  It  is  also,  however,  considered  the  low  flow  period. 
This  apparent  contradiction  is  a  result  of  evaporation  rates  exceeding  rainfall 
rates.  To  account  for  this  gain  and  loss,  it  is  necessary  to  compute  the  net  flow 
into  the  system  due  to  rain  and  evaporation. 

Flow  in  and  out  of  the  mesh  caused  by  rainfall/evaporation  is  accomplished 
by  specifying  a  rate  (m/sec)  into  (positive)  or  out  of  (negative)  the  specific 
elements.  The  inflow  caused  by  rainfall  was  computed  converting  the  recorded 
rainfall  rate  in  units  of  inches/day  into  meters/second.  Rainfall  intensity  measured 
at  thirteen  District  sites  in  the  St.  Johns  River  Basin  varied  from  site  to  site,  but 
the  frequency  of  rainfall  events  for  all  stations  and  their  relative  intensity  within 
each  station  were  similar.  The  same  fact  was  observed  for  wind  speed  and 
direction  measured  at  six  sites.  These  observations  demonstrate  that  the  major 
meteorological  events  affecting  the  behavior  of  the  Lower  St.  Johns  River  and 
Estuary  are  regional  and  can  be  assumed  uniform  over  the  entire  model.  Daily 
median  values  of  rainfall  from  among  the  13  District  sites  were  used.  Hourly 
wind  data  used  were  measured  at  Shands  Bridge,  with  missing  records  filled  using 
Dames  Point  data.  Shands  Bridge  was  selected,  since  that  is  where  the  river  is  at 
its  widest.  Monthly  evaporation  rates  were  obtained  from  a  NOAA  gauge  located 
in  Gainesville.  A  pan  correction  factor  of  0.78  was  used. 

Figures  12  and  13  show  the  rainfall  and  evaporation  data  for  the  time  period 
under  study.  Volumes  were  obtained  by  multiplying  the  flow  rate  by  the  area  of 
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the  river,  excluding  the  ocean  and  1-D  elements.  Figure  14  illustrates  the 
cumulative  flow  into  and  out  of  the  system  during  the  study  period.  The  low  flow 
and  high  flow  periods  are  easily  delineated.  Wind  speed  and  direction  are  shown 
in  Figure  15. 


Marsh  Approximation  with  One-Dimensional 
Elements 


To  increase  computer  efficiency,  some  marsh  areas  and  streams  were 
approximated  using  one-dimensional  elements.  These  areas  included: 

a.  Pablo  Creek  and  marshes, 

b.  Chicopit  Bay, 

c.  Greenfield  Creek, 

d.  Sisters  Creek  and  marshes 

e.  Cedar  Point  Creek  and  marshes, 

f.  Ribault  River, 

g.  Moncrief  Creek, 

h.  Parts  of  the  Ortega  River  and  Fishing  Creek, 

i.  Part  of  Arlington  River,  Pottsburg  Creek,  Little  Pottsburg 
Creek,  and  Strawberry  Creek, 

j.  Black  Creek, 

k.  Murphy  Creek, 

l.  Dunns  Creek, 

m.  St.  Johns  River  between  Lake  George  and  Buffalo  Fluff,  and 

n.  Lake  George 
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Field  Data 


Tide  and  salinity  data  for  verification  were  supplied  by  SJRWMD.  These 
data  cover  a  period  of  approximately  300  days,  starting  1  April  1996.  The  data  set 
included  hourly  tidal  elevation  data  at  nine  locations  and  hourly  temperature  and 
conductivity  (used  to  compute  salinity)  data  at  five  locations. 

The  nine  tidal  stations  were  located  at  the  following  locations,  which  span  the 
length  of  the  LSJRE  and  provide  a  continuous  set  of  verification  data  (see  Figure 
16):  Mayport,  Dames  Point,  U.S.  Army  Corps  of  Engineer  Dredge  Depot,  Main 
Street  Bridge,  Buckman  Bridge,  Shands  Bridge,  Racy  Point,  Palatka,  and  Buffalo 
Bluff.  The  salinity  stations  were  located  at  Dames  Point,  Acosta  Bridge, 

Buckman  Bridge,  Shands  Bridge,  and  Buffalo  Bluff  (see  Figure  17). 

The  salinity  data  were  collected  at  the  surface,  mid-depth,  and  bottom  for  all 
stations  but  Buffalo  Bluff  where  only  surface  data  were  collected.  Surface  and 
bottom  data  were  used  for  the  verification.  Gaps  in  the  surface  salinity  data  were: 
Buffalo  Bluff  from  1  July  to  8  August,  Dames  Point  from  24  June  to  2  July,  and 
Acosta  Bridge  from  3  July  to  27  August.  Gaps  in  the  bottom  salinity  data  were: 
Shands  Bridge  from  29  May  to  29  August,  and  Dames  Point  from  13  July  to  28 
August.  Measurements  at  Shands  Bridge  and  Buffalo  Bluff  represent  salinity 
values  less  than  1  ppt  having  no  identifying  shape,  form,  or  structure. 

The  model  was  verified  to  water  surface  elevation  and  water  discharge  by 
adjusting  the  mesh  and  by  selecting  the  model  parameters  of  Manning's  roughness 
coefficient  (n)  and  horizontal  eddy  viscosity.  To  provide  adequate  flows  for 
verification,  discharge  data  were  provided  by  the  SJRWMD  for  eight  locations 
(see  Figure  18).  These  eight  locations  are  located  near  or  at:  Bayard  Point,  Corps 
Dock,  Dames  Point,  Doctors  Inlet,  Main  Street,  Mayport  (Marker  22),  Blount 
Island  Channel  (Marker  36),  and  Fulton-Dame  Cutoff  Range  (Marker  38). 
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3-D  Model  Verification 


The  model  was  operated  in  three-dimensional  mode.  The  three-dimensional 
mode  uses  one-,  two-,  and  three-dimensional  elements. 

The  3-D  model  was  run  for  a  total  of  150  days.  The  first  30  days  represented 
the  period  29  April  to  28  May  1996  and  served  as  the  spin-up  period.  Spin-up  is 
necessary,  because  the  initial  conditions  are  not  precisely  known.  In  studies 
involving  hydrodynamics  only,  the  spin-up  period  can  consist  of  as  little  as  one  or 
two  days.  With  hydrodynamics  and  transport  (salinity),  however,  it  can  require 
weeks  to  months  of  simulation  to  reach  equilibrium,  real-world  conditions  of 
salinity.  The  initial  salinity  field  for  the  3-D  simulation  was  approximated  by 
interpolating  average  field  salinity  values  to  the  numerical  model  mesh.  These 
values  were  input  by  RMA10-WES  and  the  salinity  values  assigned  to  the 
corresponding  3-D  numerical  mesh  node.  From  that  starting  point,  RMA10-WES 
computed  the  hydrodynamics  and  salinities. 


Tidal  elevation 

Figure  16  shows  the  location  of  the  nine  prototype  tidal  stations.  Water 
surface  elevation  data  collected  at  Buckman  Bridge,  Dames  Point,  USCOE 
Dredge  Depot,  Palatka,  and  Shands  Bridge  contained  various  gaps.  Figures  17-25 
show  the  prototype  data  plotted  with  the  3-D  model  results  for  the  nine  stations 
(all  locations  have  been  plotted  with  reference  to  NGVD). 

Very  good  agreement  was  observed  at  the  Mayport  station  for  tidal  range  and 
phasing.  Moving  upstream,  the  water  surface  elevation  computed  by  the  model 
tends  to  remain  higher  than  the  prototype.  Several  unsuccessful  efforts  to  lower 
the  water  surface  elevation  of  the  model  took  place  by  modifying  numerous  model 
parameters.  Since  the  model  salinity  was  observed  to  be  higher  than  the  prototype 
downstream  from  Jacksonville  downtown,  an  investigation  of  the  model  was 
performed  by  uncoupling  hydrodynamic  calculations  to  the  salinity.  A  good 
match  of  the  water  surface  elevation  between  model  and  prototype  was  obtained 
as  a  result  of  this  change  (see  Figure  26).  The  explanation  for  this  model  behavior 
is  that  the  density  of  the  fresh  water  located  in  the  upstream  part  of  the  model  is 
lower  than  the  density  of  the  salty  water  in  the  downstream  section  of  the  LSJRE 
model.  Therefore,  the  water  surface  elevation  upstream  is  forced  to  increase  to 
mass  balance  against  heavier  water  located  downstream.  The  tidal  range 
throughout  the  model  was  also  improved  by  uncoupling  salinity.  After  these 
findings,  verification  efforts  were  concentrated  on  improving  salinity  results. 

The  water  level  datum  planes  add  uncertainty  to  the  verification.  They  are 
notoriously  unreliable  and  could  be  responsible  for  some  of  the  difference 
between  prototype  data  and  model  results. 
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Water  flow 


Figure  43  shows  the  location  of  eight  discharge  measurement  sites  as  well  as 
the  dates  when  the  measurements  took  place.  Figures  44-5 1  show  the  prototype 
discharge  plotted  with  the  3-D  model  results. 

The  results  at  Bayard  Point  show  good  agreement  in  both  magnitude  and 
phasing.  Although  the  magnitude  of  the  water  flowing  through  the  Blount  Island 
Channel  (Marker  36)  is  well  matched  by  the  model,  there  is  an  hour  lag  in  the 
phase.  The  difference  in  phase  is  due  to  the  geometry  simplification  of  the 
adjacent  marshes  and  the  lack  of  data  to  verify  the  timing  of  the  water  flow 
interchange  between  the  marshes  and  the  main  stem  of  the  river.  The  magnitude 
of  the  discharge  at  Mayport  could  have  also  been  overestimated  by  the 
simplification  of  the  marshes.  Another  possible  factor  overestimating  the 
discharge  magnitude  at  Mayport  is  the  simplified  geometry  of  the  river  entrance 
where  jetties  were  substituted  by  a  smooth  rounded  entrance. 

The  model  and  prototype  discharge  data  at  the  other  five  locations  match  well. 
Consistently,  the  magnitude  of  the  model  discharges  was  slightly  higher.  This 
observation  agrees  with  the  fact  that  water  surface  elevations  were  also  higher 
during  the  simulation  period.  Figure  52  shows  how  discharges  are  affected  by 
uncoupling  salinity  to  hydrodynamic  computations. 

Attempts  to  reduce  discharges  by  increasing  the  bed  roughness  were 
unsatisfactory  because  the  tide  range  throughout  the  model  was  compromised  as  a 
result. 


Salinity 

Figure  27  shows  the  location  of  the  salinity  stations.  Data  were  collected  at 
three  depths  at  all  of  the  stations  except  Buffalo  Bluff.  The  surface  and  bottom 
salinity  values  for  model  and  prototype  are  shown  for  four  of  the  stations  in 
Figures  28-40.  Buffalo  Bluff  only  shows  surface  salinity. 

Referring  to  the  Dames  Point,  Acosta  Bridge,  and  Buckman  Bridge  data,  the 
three-dimensional  model  captures  the  salinity  trend,  including  isolated  peaks.  A 
good  example  is  observed  on  Julian  day  194.  All  salinity  stations  show  a  water 
surface  elevation  and  salinity  peak  in  prototype  data  and  model  results.  Figure  14 
shows  an  increase  in  wind  speed  as  well  as  a  change  in  direction  on  the  same  day. 
Winds  were  blowing  northeast  and  suddenly  changed  to  southwest.  The  response 
of  the  model  justifies  the  assumption  of  using  uniform  winds  for  this  study.  At 
Shands  Bridge  and  Buffalo  Bluff,  the  model  results  represent  well  the  prototype 
conditions  with  salinity  values  less  than  1  ppt  and  with  no  identifying  shape,  form, 
or  structure. 

As  mentioned  earlier,  model  salinities  are  somewhat  higher  than  those 
observed  in  the  lower  (seaward)  half  of  the  model.  Two  aspects  of  the  model  are 
considered  as  likely  reasons  for  the  higher  average  salinities  computed  in  the 


lower  portion  of  the  model:  a)  ungauged  watershed  runoff  and  groundwater  flow 
into  the  river;  and  b)  simplified  geometry,  particularly  at  the  river  entrance.  Since 
both  model  aspects  have  an  effect  on  the  river  freshwater  levels,  additional 
simulations  were  performed  to  test  the  sensitivity  of  the  model  to  higher 
freshwater  inflows.  Results  of  one  of  the  simulations  are  shown  in  Figure  41 .  In 
this  simulation,  130.26  m3/sec  (4,600  cfs)  of  fresh  water  were  distributed  over  the 
surface  of  the  river.  This  amount  of  water  yields  a  satisfactory  fit  between 
prototype  and  modeled  salinity  values  at  Acosta  Bridge.  Although  the  application 
of  the  additional  fresh  water  dropped  salinity  considerably  over  the  entire  model,  a 
different  distribution  and  magnitude  of  freshwater  inflow  could  yield  better 
results.  Again,  the  simplification  in  the  geometry  of  the  St.  Johns  River  entrance 
could  overestimate  the  interchange  of  water  in  the  downstream  section  of  the 
river.  Also,  the  secondary,  northern  entrance  was  not  modeled,  and  the 
surrounding  marshes  were  simplified  in  the  model. 

The  model  still  showed  the  same  stratification  trends  found  in  the  prototype, 
although  more  markedly.  Stratification  was  found  in  the  model  as  far  upstream  as 
Buckman  Bridge,  reproducing  the  same  prototype  behavior. 

We  interpret  these  salinity  results  to  indicate  that  the  model  successfully 
reproduces  the  overall  salinity  trends  within  the  LSJRE,  with  the  understanding 
that  some  local  areas  may  not  be  as  well  reproduced  as  the  five  stations  shown 
until  verification  of  those  areas  is  checked  and  any  necessary  adjustments  made. 
The  model  simulates  salinity  pulses  caused  by  changes  in  freshwater  flow  and 
stratification  changes,  but  the  tidal  surface  elevation  throughout  the  model  is  the 
variable  playing  a  major  roll  in  the  overall  salinity  behavior,  as  shown  in  Figure 
42. 


Computational  Speed 


The  3-D  geometry  used  to  generate  the  results  discussed  is  computationally 
time  intensive,  but  the  model  solution  converges  with  the  use  of  only  three 
iterations  per  time  step,  two  for  hydrodynamic  calculations  and  one  for  salinity. 
Half-hour  time  steps  were  used  for  this  study. 

The  platform  used  to  run  the  model  is  an  SGI  Origin  2000.  This  is  part  of  the 
ERDC  High  Performance  Computing  Center.  The  SGI  Origin  is  a  cc-NUMA 
distributed  shared  memory  architecture  that  contains  1 12  processors,  64  Gbytes  of 
memory,  and  500  Gbytes  of  hard  disk  storage.  We  simulated  as  much  as  10  days 
during  12  wall  clock  hours,  or  28  sec  per  iteration. 
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4  Discussion  and  Conclusions 


Discussion 


The  purpose  of  a  numerical  model  is  to  achieve  insight  and  understanding 
about  the  system  being  modeled  to  help  in  a  decision  making  process.  A  complete 
understanding  of  a  system  cannot  be  obtained  within  ordinary  constraints  of  time 
and  money;  therefore,  a  model  must  be  designed  to  satisfy  specific  criteria  in 
anticipation  of  its  intended  uses.  For  example,  if  a  LSJRE  model  will  be  used  to 
study  the  salinity  and  sedimentation  effects  of  a  navigation  channel  enlargement, 
the  modeler  will  build  a  different  model  than  if  its  purpose  is  to  study  the  effects 
of  non-point  pollution  on  water  quality.  The  former  will  portray  very  high 
resolution  in  the  navigation  channel,  while  the  latter  may  ignore  the  exact  location 
of  the  channel  as  long  as  the  river  cross  section  is  adequately  described. 

The  two  aspects  of  waterway  modeling  most  important  for  faithful 
reproduction  of  hydrodynamic  and  transport  processes  are  an  accurate  rendering 
of  the  geometry  (bank  lines  and  bathymetry)  and  accurate  specification  of  the 
boundary  conditions  (tides,  river  flows,  etc.). 

The  model  described  here  was  designed  to  provide  a  satisfactory  reproduction 
of  estuarine  hydraulics  and  salinity  transport.  This  objective  was  accomplished. 
With  presently  available  micro  and  mini  computers,  the  model  runs  too  slowly  to 
be  considered  practical  for  continuous  simulation  periods  longer  than  a  few 
months.  The  use  of  a  mainframe  supercomputer  is  necessary  for  the  feasibility  of 
a  long-term  simulation  study.  Lower  resolution  versions  with  a  more  structured 
layout  (i.e.,  a  nearly  uniform  number  of  elements  across  the  estuary)  will  run 
considerably  faster.  For  water  quality  modeling,  which  typically  requires  multi¬ 
year  simulations,  such  a  lower  resolution  and,  thus,  a  computationally  more 
efficient  version  of  this  model  is  adequate  and  should  be  considered. 

Accurate  specification  of  boundary  conditions  is  often  more  difficult  than 
accurate  specification  of  system  geometry.  Water  level  datum  planes  are 
notoriously  unreliable,  and  they  should  always  be  subjected  to  critical  evaluation. 
Even  routine  discharge  measurements  are  less  precise  than  desired.  In  this  model, 
we  added  rain  and  evaporation  estimates  to  achieve  salinity  verification,  but  did 
not  include  groundwater  flows  or  local  runoff,  which  can  significantly  affect 
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salinity  and  mixing. 


Future  use  of  the  model  should  bear  these  limitations  in  mind  and  (a)  revise 
the  mesh  if  new  uses  require  it,  (b)  recheck  the  local  verification  with  high  quality 
data  if  specific  zones  are  of  high  interest,  and  (c)  refine  boundary  conditions  of  the 
model,  depending  on  the  usage. 


Conclusions 


A  hydrodynamics  and  salinity  transport  model  of  the  Lower  St.  Johns 
River  and  Estuary  has  been  constructed  and  verified.  The  RM A 10- WES 
model  described  here  satisfactorily  reproduces  tidal  elevations,  discharges, 
and  salinity  trends. 

Although  the  model  verification  is  satisfactory,  improvements  can  be  made: 

a.  More  comprehensive,  synoptic  field  data  collected  specifically  for  the 
purpose  of  modeling  is  needed  to  check  local  verification  if  detailed 
results  are  desired  at  locations  among  the  verification  stations  used 
herein. 

b.  Improved  definition  of  boundary  conditions,  including  local  runoff  as 
well  as  groundwater  interactions  will  improve  the  model. 

c.  A  mesh  geometry  more  closely  representing  the  river  entrance  and  the 
extensive  marshes  located  to  the  north  of  the  LSJRE  is  necessary. 
Prototype  data  (water  velocity  and  discharge)  collected  within  the 
marshes  and  their  connection  to  the  main  stem  of  the  river  should 
complement  improvements  in  the  geometry. 

The  model  is  suitable  for  addressing  questions  of  general  hydrodynamics  and 
salinity  intrusion  within  the  main  body  of  the  LSJRE.  Off-channel  tributaries  and 
marshes  have  been  modeled  for  their  effect  on  the  main  waterway,  but  have  not 
been  verified  to  reproduce  local,  internal  hydrodynamic  behavior.  For  future 
applications,  the  level  of  verification,  degree  of  resolution  and  appropriateness  of 
boundary  conditions  should  be  carefully  evaluated  in  each  case  to  ensure  a  proper 
match  between  uses  and  model  behavior  in  the  areas  of  interest. 

The  model  is  computationally  intensive.  Simulations  longer  than  a  few 
months  can  be  achieved  with  the  use  of  a  mainframe  supercomputer  or  by 
optimizing  the  mesh  for  balanced  accuracy  and  speed,  matching  mesh  resolution 
to  the  questions  posed. 
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Table  1 

Turbulent  Exchange  Coefficients  and 


Material 


4788  (100.00) 


11970  (250.00) 


9576  (200.00) 


9576  (200.00) 


11970(250.00) 


4788(100.00) 


4788  (100.00) 


3591  (75.00) 


11970(250.00) 


9576  (200.00) 


11970(250.00) 


11970(250.00) 


11970  (250.00) 


3591  (75.00) 


2394  (50.00) 


1197(25.00) 


3591  (75.00) 


4788  (100.00) 


47880  (1000.00) 


35910  (750.00) 


14364  (300.00) 


2394  (50.00) 


1197  (25.00) 


|  2394  (50.00) 


2394  (50.00) 


2394  (50.00) 


2394  (50.00) 


2394  (50.00) 


2394  (50.00) 


2394  (50.00) 


2394  (50.00) 


2394  (50.00) 


|  2394  (50.00) 


2394  (50.00) 


2394  (50.00) 


2394  (50.00) 


2394  (50.00) 


2394  (50.00) 


1 2394  (50.00) 


3591  (75.00) 


3591  (75.00) 


3591  (75.00) 


3591  (75.00) 


3591  (75.00) 


I  3591  (75.00) 


3591  (75.00) 


3591  (75.00) 


3591  (75.00) 


3591  (75.00) 


3591  (75.00) 


3591  (75.00) 


3591  (75.00) 


3591  (75.00) 


3591  (75.00) 


1 3591  (75.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


|  5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


1 5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


5.57  (60.00) 


Manning’s  n  Values 


Vertical  Manning’s 

Diffusion  n  Value 

m2/s  (ft2/s) 


0.0000929(0.001)  0.018 


0.00929  (0.100) 


0.0000929(0.001)  0.018 


0.0000929(0.001)  0.018 


0.0000929(0.001)  0.018 


0.0000929(0.001)  0.018 


0.0000929(0.001)  0.024 


0.0000929  (0.001 )  0.028 


0.0000929(0.001)  0.022 


0.0000929(0.001)  1 0.018 


0.00929  (0.100) 


0.0000929(0.001)  10.024 


0.0000929(0.001)  0.018 


0.0000929(0.001)  0.018 


0.0000929(0.001)  0.018 


0.0000929  (0.001)  0.026 


0.0000929(0.001)  0.012 


0.0000929(0.001)  1 0.024 


0.01858  (0.200) 


0.01858  (0.200) 


0.013935  (0.150) 


Number  of 
3D  Layers 


0.024 


0.024 


0.024 


0.018 


0.0000929  (0.001 )  0.026 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.024 


0.024 


0.0000929(0.001)  0.024 


0.0000929(0.001)  0.012 


0.0000929(0.001)  0.012 


0.0000929(0.001)  0.012 


0.012 


0.0000929(0.001)  0.012 


0.0000929(0.001)  0.012 


0.0000929(0.001)  0.012 


0.0000929(0.001)  0.012 


0.0000929(0.001)  0.012 


0.0000929(0.001)  0.012 


0.0000929(0.001)  0.012 


0.0000929(0.001)  0.012 


0.0000929  (0.001) 


0.0000929(0.001)  0.012 


0.0000929(0.001)  0.012 


0.012 


0.0000929(0.001)  0.012 


0.0000929(0.001)  0.012 


a.  Mesh  boundary  and  boundary  condition  locations  (Sheet  1  of  5) 
Figure  1 .  Lower  St.  Johns  River  and  Estuary  Model  mesh 


b.  Lower  portion  of  mesh 
Figure  1 .  (Sheet  2  of  5) 


c.  Lower  middle  portion  of  mesh 
Figure  1 .  (Sheet  3  of  5) 


d.  Upper  middle  portion  of  mesh 
Figure  1 .  (Sheet  4  of  5) 


e.  Upper  portion  of  mesh 


Figure  1.  (Sheet  5  of  5) 


Figure  2.  Cumulative  inflow  from  tributaries  into  the  Lower  St.  Johns  River  and 
Estuary  (To  convert  volume  to  cubic  meters,  multiply  by  0.02832) 
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Figure  3.  Big  Davis  Creek  daily  average  flow  (To  convert  discharge  to  cubic 
meters  per  second,  multiply  by  0.02832) 


DISCHARGE,  CFS 


Figure  4.  Black  Creek  daily  average  flow  (To  convert  discharge  to  cubic  meters  per 
second,  multiply  by  0.02832) 


Figure  5.  Buffalo  Bluff  daily  average  flow  (To  convert  discharge  to  cubic  meters 
per  second,  multiply  by  0.02832) 


Figure  6.  Deep  Creek  daily  average  flow  (To  convert  discharge  to  cubic  meters  per 
second,  multiply  by  0.02832) 


Figure  7.  Dunns  Creek  daily  average  flow  (To  convert  discharge  to  cubic  meters 
per  second,  multiply  by  0.02832) 


Figure  8.  Ortega  River  daily  average  flow  (To  convert  discharge  to  cubic  meters 
per  second,  multiply  by  0.02832) 


Figure  9.  Rice  Creek  daily  average  flow  (To  convert  discharge  to  cubic  meters  per 
second,  multiply  by  0.02832) 


a.  Daily  average  flow 


b.  Corrected  daily  average  flow 

Figure  10.  Cedar  River  average  flow  (To  convert  discharge  to  cubic  meters  per 
second,  multiply  by  0.02832) 
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b.  Corrected  daily  average  flow 


Figure  11.  Six  Mile  Creek  average  flow  (To  convert  discharge  to  cubic  meters  per 
second,  multiply  by  0.02832) 


EVAPORATION,  IN 


Figure  12.  Rainfall  (median  values  among  13  district  sites)  (To  convert  rainfall  to 
mm/day,  multiply  by  25.4) 


Figure  13.  Evaporation  at  Gainesville,  1996  (To  convert  evaporation  to  mm, 
multiply  by  25.4) 
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14.  Cumulative  rainfall/evaporation  into  the  Lower  St.  Johns  River  and 
Estuaiy  (To  convert  volume  to  m3,  multiply  by  0.0283 1 685) 


Figure  15.  Wind  (Shands  Bridge,  missing  records  filled  using  Dames  Point) 
(To  convert  magnitude  to  km,  multiply  by  1.609) 


Figure  16.  LSJRE  tide  stations 
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Figure  17.  Water  surface  elevation  at  Mayport  (Continued)  (To  convert  elevation  to 
meters,  multiply  by  0.3048) 
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Figure  17.  (Concluded) 
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Figure  18.  Water  surface  elevation  at  Dames  Point  (Continued)  (To  convert  elevation 
to  meters,  multiply  by  0.3048) 
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Figure  18.  (Concluded) 
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Figure  19.  Water  surface  elevation  at  USCOE  Dredge  Depot  (Continued)  (To  convert 
elevation  to  meters,  multiply  by  0.3048) 
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Figure  19.  (Concluded) 
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Figure  20.  Water  surface  elevation  at  Main  Street  Bridge  (Continued)  (To  convert  elevation 
to  meters,  multiply  by  0.3048) 
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Figure  20.  (Concluded) 
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Figure  21.  Water  surface  elevation  at  Buckman  Bridge  (Continued)  (To  convert  elevation  to 
meters,  multiply  by  0.3048) 
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Figure  21.  (Concluded) 
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Figure  22.  Water  surface  elevation  at  Shands  Bridge  (Continued)  (To  convert  elevation  to 
meters,  multiply  by  0.3048) 
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Figure  22.  (Concluded) 
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Figure  23.  Water  surface  elevation  at  Racy  Point  (Continued)  (To  convert  elevation  to 
meters,  multiply  by  0.3048) 
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Figure  23.  (Concluded) 
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Figure  24.  Water  surface  elevation  at  Palatka  (Continued)  (To  convert  elevation  to  meters, 
multiply  by  0.3048) 
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Figure  24.  (Concluded) 
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Figure  25.  Water  surface  elevation  at  Buffalo  Bluff  (Continued)  (To  convert  elevation  to 
meters,  multiply  by  0.3048) 
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Figure  25.  (Concluded) 
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Figure  26.  Effect  of  uncoupling  hydrodynamic  and  salinity  computations  on  water  surface 
elevation  at  Shands  Bridge  (To  convert  elevation  to  meters,  multiply  by  0.3048) 
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Figure  27.  LSJRE  salinity  stations 


SALINITY,  PPT 


JULIAN  DAY,  1996 


PROTOTYPE  SURFACE  . PROTOTYPE  BOTTOM 


Figure  28.  Prototype  salinity  at  Dames  Point  (Continued) 
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Figure  28.  (Concluded) 
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Figure  29.  Model  salinity  at  Dames  Point  (Continued) 
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Figure  29.  (Concluded) 
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Figure  30.  Surface  salinity  at  Dames  Point  (Continued) 


Figure  30.  (Concluded) 
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Figure  31 .  Bottom  salinity  at  Dames  Point  (Continued) 
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Figure  31.  (Concluded) 
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Figure  32.  Prototype  salinity  at  Acosta  Bridge  (Continued) 
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Figure  32.  (Concluded) 
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Figure  33.  Model  salinity  at  Acosta  Bridge  (Continued) 
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Figure  33.  (Concluded) 


Figure  34.  Surface  salinity  at  Acosta  Bridge  (Continued) 


Figure  34.  (Concluded) 
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Figure  35.  Bottom  Salinity  at  Acosta  Bridge  (Continued) 
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Figure  35.  (Concluded) 
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Figure  36.  Prototype  salinity  at  Buckman  Bridge  (Continued) 
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Figure  36.  (Concluded) 
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Figure  37.  Model  salinity  at  Buckman  Bridge  (Continued) 
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Figure  37.  (Concluded) 
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Figure  38.  Surface  salinity  at  Buckman  Bridge  (Continued) 
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Figure  38.  (Concluded) 
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Figure  39.  Bottom  salinity  at  Buckman  Bridge  (Continued) 
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Figure  39.  (Concluded) 
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a.  Shands  Bridge 
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b.  Buffalo  Bluff 


Figure  40.  Salinity  at  Shands  Bridge  and  Buffalo  Bluff 
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Figure  42.  Relationship  between  daily  averages  of  salinity  and  water  surface  elevation 
(Acosta  Bridge) 
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Figure  43.  Location  of  water  discharge  measurements 


Figure  44.  Water  discharge  at  Mayport  (Marker  22)  (To  convert  dicharge  to  cubic  meters 
per  second,  multiply  by  0.02832) 
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Figure  45.  Water  discharge  at  Fulton-Dame  Point  Cutoff  Range  (Marker  38)  (To  convert 
discharge  to  cubic  meters  per  second,  multiply  by  0.02832) 
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46.  Water  discharge  at  Blount  Island  Channel  (Marker  36)  (To  convert  discharge 
to  cubic  meters  per  second,  multiply  by  0.02832) 
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Figure  47.  Water  discharge  at  Dames  Point  (To  convert  discharge  to  cubic  meters  per 
second,  multiply  by  0.02832) 


Figure  48.  Water  discharge  at  Corps  Dock  (To  convert  discharge  to  cubic  meters  per 
second,  multiply  by  0.02832) 


Figure  49.  Water  discharge  at  Main  Street  (To  convert  discharge  to  cubic  meters  per 
second,  multiply  by  0.02832) 
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Figure  50.  Water  discharge  at  Doctors  Inlet  {To  convert  discharge  to  cubic  meters  per 
second,  multiply  by  0.02832) 


Figure  51 .  Water  discharge  at  Bayard  Point  (To  convert  discharge  to  cubic  meters  per 
second,  multiply  by  0.02832) 
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Figure  52.  Effect  of  uncoupling  hydrodynamic  and  salinity  computations  on  water  discharge 
at  Doctors  Inlet  (To  convert  discharge  to  cubic  meters  per  second,  multiply  by 
0.02832) 


Appendix  A 

The  Finite  Elements  Approach 
to  Modeling  Hydrodynamics 
and  Sediment  Transport  Using 
RMA10-WES 


Introduction 

A  model  of  a  particular  estuary,  river  or  reservoir,  is  composed  of  several  parts. 
There  is  the  geometric  description  which  not  only  includes  the  X,  Y,  and  Z 
locations  of  points  chosen,  but  also  the  delineation  of  bed  cover  such  as  grasses, 
sands,  etc.  There  is  also  the  specification  of  boundary  conditions  including  all  of 
the  forcing  mechanisms,  tides,  wind,  inflow,  and  any  constituents.  The  third  part 
is  numerical  code.  This  code  has  at  its  core,  the  basic  laws  of  physics,  in 
particular  the  relationship  between  force  and  acceleration.  A  consistent  model 
code  is  a  faithful  representation  of  these  equations,  for  sufficient  resolution.  Since 
these  equations  are  universal,  the  range  of  applicability  of  the  code  is  limited  by 
the  additional  empirical  relationships  to  describe  turbulent  processes,  bed 
roughness,  equations  of  state,  and  wind  stress,  and  by  any  simplifying 
assumptions  made  in  the  equations. 

There  is  some  degree  of  uncertainty  associated  with  the  geometric  description, 
the  boundary  conditions  and  in  the  numerical  code.  Often  the  error  indicated  by 
these  uncertainties  is  highest  for  the  specification  of  the  boundary  conditions, 
followed  by  the  geometric  description,  with  the  numerical  code  providing  the  least 
of  the  error.  The  precise  values  of  inflow  and  wind  fields,  particularly  in  estuarine 
models,  are  usually  fairly  crude.  The  inflows,  for  example,  are  determined  by 
gauges  far  upstream  from  the  estuary.  They  leave  out  large  ungauged  downstream 
areas.  Also,  at  the  gauge  the  relationship  between  gauge  elevation  and  discharge 
is  far  from  perfect.  This  relationship  changes  over  time  and  flow  conditions. 

These  errors  can  be  minimized  by  careful  model  verification.  Also,  modelers 
can  use  sound  testing  procedures  to  eliminate  much  of  the  error  from  the  decision 
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process.  These  techniques  include  plans  to  base  comparisons  and  sensitivity 
analysis. 

In  this  paper,  one  hydrodynamic  code  is  described,  RMA10-WES.  The 
Galerkin-based  finite  element  model  RMA10-WES  is  an  ERDC  adaptation  of  the 
RMA-10  code  developed  by  King  (1993)'.  This  code  computes  time-varying 
open-  channel  flow  and  salinity/temperature  transport  in  one,  two,  and  three 
dimensions.  It  invokes  the  hydrostatic  pressure  and  mild  slope  assumption. 
Vertical  turbulence  is  supplied  using  a  Mellor-Yamada  Level  II  (Mellor  and 
Yamada  1982)  k-(  approach  modified  for  stratification  by  the  method  of 
Henderson-Sellers  (1984).  The  salinity/density  relationship  is  based  upon 
Pritchard  (1982).  The  RMA2  code  is  a  predecessor  to  RMA10-WES  but  is 
limited  to  at  most  two  dimensional  cases.  The  formulation  is  similar  to  RMA10- 
WES  and  so  the  only  descriptions  reported  herein  will  be  those  of  RMA10-WES. 


Model  Description 


The  full  three-dimensional  equations  are  reduced  to  a  set  of  two  momentum 
equations,  an  integrated  continuity  equation,  a  convection-diffusion  equation,  and 
an  equation  of  state.  The  simplification  is  a  result  of  the  hydrostatic  pressure 
approximation: 


Du 
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ds 
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(3) 


(4) 


P  =  F(s) 


(5) 


Elevation-related  terms  are  defined  in  Figure  1. 


1  References  cited  in  this  appendix  are  listed  in  the  References  at  the  end  of  the  main  text. 
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where 


1 

' 

Exx  du 

Eyx  5V 

dx 

dx 

Exy  5u 

’  ;  <J  y  -• 

Eyy  5V 

dy 

dy 

Ex2 

Eyz  5V 

dz  . 

dz  . 

and 

p  =  density 

u,v,w  =  x,y,z  velocity  components 
t  =  time 
P  =  pressure 

rx  =  pQ.v-^Ua(ua  Va)  +ipw2  cos(G) 
C 

ry  =  -PaU-PSVa(uai  Va)  +y/w2s\n( Q) 
C 

Q  =  2cosin(^) 

co  =  rate  of  angular  rotation  of  the  earth 
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<j>  =  local  latitude 
g  =  gravitational  acceleration 
C  =  Chezy  or  Manning  friction  formulation 
if/  =  a  coefficient  from  Wu  (1980) 

W  =  wind  speed 

0  =  wind  direction  counterclockwise  from  easterly 
h  =  depth 

uq,\>q  -  x,y  velocity  components  at  the  water  surface 
C  =  water  surface  elevation 
ua,va  =  x,y  velocity 
a  =  bed  elevation 
s  =  salinity 

DxJDy,D .  =  diffusion  coefficient  for  salt 
E  =  eddy  viscosity  components 


The  continuity  equation: 


du  dv 
—  +  —  + 

dx  dy 


(6) 


is  solved  as  a  second  part  of  each  solution  step.  Equation  6  is  converted  to  an 
appropriate  boundary  value  problem  through  differentiation  with  respect  to  z. 
After  rearrangement  it  takes  the  form 


d2w  _  d  f  du  dv N 
dz2  dz\dx  dy  y 


(7) 


subject  to  boundary  conditions  specified  for  the  water  surface  and  the  bed: 


W£-U£ 


dx 


Vf 


d£_+dh 
dy  dt 


at  the  water  suface 


(8) 


and 
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at  the  bed . 


(9) 


Wa  ~Ua 


da 

dy 


Note  that  in  these  equations  the  values  of  u  and  v  will  be  known  at  all  locations 
from  the  previous  part  of  the  solution  step.  Values  of  w  in  this  solution  are  used 
in  the  next  iteration  for  u,  v,  h,  and  s. 

The  geometric  system  varies  with  time;  i.e.,  the  water  depth  h  varies  during  the 
simulation.  To  develop  an  Eulerian  form  for  the  solution,  it  is  desirable  to 
transform  this  system  to  one  that  can  be  described  with  a  constant  geometric 
structure.  Early  development  of  the  model  (King  1982)  used  a  cr-transformation 
in  which  the  bed  and  the  water  surface  are  transformed  to  constants.  In  a  later 
analysis  of  this  method,  King  (1985)  pointed  out  that  at  locations  where  a  sharp 
break  in  bottom  profile  occurs,  the  transformation  is  not  unique  and  momentum  in 
the  component  directions  may  not  be  correctly  preserved.  An  alternative 
transformation  that  preserves  the  bottom  profile  as  defined,  but  transforms  the 
water  surface  to  a  constant  elevation  is  now  used  {z7 transformation). 

This  transformation  is  defined  by: 

*-=  x  (10) 


y~=  y 


(ii) 


z  -  =  a  +  (z  -  a)  — - — —  (12) 

h 

where  b  is  the  fixed  vertical  location  to  which  the  water  surface  will  be  trans¬ 
formed.  Equations  1-6  and  7-9  then  incorporate  the  transformation  (Equations 
10-12). 

Another  advantage  of  this  transformation  is  that  it  produces  z7  =  constant  lines 
that  are  close  to  horizontal,  i.e.,  z  =  constant  lines.  This  results  in  less  fictitious 
density-driven  currents  near  bed  profile  breaks  (Stelling  and  van  Kester  1993). 

Since  stratification-related  phenomena  are  usually  nearly  horizontal,  it  is 
important  that  the  transformation  leave  constant  surfaces  horizontal.  Considering 
the  pressure  gradient  (due  to  the  density  gradient)  in  this  transformation  it 
produces: 


dP  _  dP  dP  dz- 
dx  dxr  dz-  dx 


(13) 


In  a  strongly  stratified  stagnant  system  this  pressure  gradient  should  be  zero. 
However,  note  that  Equation  13  in  the  transformed  system  is  dependent  upon  two 
terms  (each  of  which  could  be  large)  to  cancel  each  other.  This  could  cause 
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artificial  currents  due  to  truncation  and  roundoff  error.  A  transformation  in  which 
dzv!dx  ~  0,  i.e.,  z*  ~  z,  will  reduce  this  problem.-  Figure  2  shows  an  example  in 
which  a  12.2-m-deep  (40-ft-deep)  channel  passes  through  a  2.4-m-deep  (8-ft- 
deep)  bay.  Here  b  is  chosen  to  be  an  elevation  of  0  and  C,  is  0.6  m  (2  ft).  Near  the 
break  in  the  bed  profile  dz7lck  is  fairly  small,  or  zv  surfaces  are  nearly  horizontal. 
Contrast  this  with  the  a  transformation  in  Figure  3.  The  <7=  constant  surfaces  are 
far  from  horizontal  along  the  channel  side  slopes.  The  truncation  and  roundoff 
errors  tend  to  drive  fictitious  currents  that  cause  the  denser  salt  water  to  leave  the 
channel.  The  zr  transformation  results  in: 

- 'f-  =  0((-b )  (14) 

OX 

whereas  the  cr transformation  is: 


(15) 


which  is  much  larger. 

The  Galerkin  finite  element  approximation  of  Equations  1-4  and  7  uses  a 
quadratic  approximation  for  u,  v,  w,  and  5  and  linear  for  h  and  P.  The 


Figure  2.  Lines  of  constant  z  near  a  significant  grade  change 
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Figure  3.  Lines  of  constant  a  near  a  significant  grade  change 

nonlinearity  is  addressed  by  Newton-Raphson  iteration  at  each  time-step. 
Generally  the  iteration  process  is  split  into  calculation  of  Equations  1-3,  then  7, 
followed  by  4.  This  sequence  is  repeated  until  sufficient  convergence  is  reached. 


Sediment  Transport 


RMA2  is  an  ancestor  of  the  RMA10-WES  code  and  is  similar  in  approach. 
That  code  was  2-D  and  1-D.  It  did  not  directly  handle  sediment  transport. 

Instead,  it  supplied  the  currents  and  water  surface  elevations  to  drive  a  separate 
sediment  transport  code,  STUDH  or  SED2D.  These  sediment  transport  codes 
include  both  cohesive  and  non-cohesive  sediment  and  the  bed  interaction 
algorithms.  By  separating  the  hydrodynamics  and  sediment  transport  into  two 
codes,  the  sedimentation  process  could  not  directly  influence  the  hydrodynamics. 
In  the  development  of  RMA10-WES,  we  have  included  the  sediment  transport  in 
the  hydrodynamic  code.  At  present  only  cohesive  sediments  are  programmed. 

For  cohesive  sediment  the  erosion  and  deposition  is  described  by  time  rates 
that  are  dependent  upon  the  bed  shear  stress.  In  RMA10-WES  the  deposition  rate 
is  described  Krone  (1962). 
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(16) 


S  =  -  V  s  C(1  ~—)forC  <  Cc 

Td 

S  =  -Vkc5 /3(1-  —)forC  >  Cc  (17) 

Td 

where: 

x  =  bed  shear  stress 

Td=  critical  shear  stress  for  deposition 

Vs=  settling  velocity 

C  =  the  concentration  near  the  bed 

Cc=  critical  concentrarion,  300  mg/l 

Vk=  the  hindered  settling  velocity,  Vs/Cc  4/3 

The  erosion  rate  is  given  by  Partheniades  (1962). 

S  =  P(--1)  (18) 

Te 


where: 


P  =  erosion  rate  constant 

xe  =  the  critical  shear  stress  for  particle  erosion. 


The  bed  algorithim  is  maintained  as  a  set  of  fluctuation  layers,  each  with  as¬ 
signed  properties.  The  concentration  of  the  sediment  is  approximated 
quadratically  over  each  element  in  the  grid. 
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